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I. INTRODUCTION 

Materials composed of hard cohesionless grains, such 
as dry sand, exhibit many remarkable properties ranging 
from cluster formation in gaseous phases, (lj,J|2|,j^] to un- 
expected flows and jets in fluid-like phases, and to 
complex organization of inhomogeneous stresses in dense 
fluid or sohd phases. 0, 0, [^] In this paper we focus on 
one feature of the solid phase that has received much 
attention: the distribution P{f) of contact forces be- 
tween grains in a system that is supporting a macroscopic 
compression or shear force. In a number of experiments 
and numerical simulations involving non-cohesive grains, 
P(/) appears to decay approximately exponentially at 
large /, 0,0, 0,0,^5 which runs counter to the naive 
expectation of a Gaussian distribution about some av- 
erage /. Several theoretical analyses of model systems 
have indicated possible explanations for an exponential 
tail,0, J^, XL JJ.I 2Mj Jfll including an analytic cal- 
culation for a special case of isostatic packings of friction- 
less disks. (2l| Still, a fundamental understanding of the 
phenomenon has not been achieved. 

To motivate the problem considered in the present 
work, we recall that for generic packings of frictional 
and/or nonspherical hard particles geometrical con- 
straints permit the formation of more contacts than 
would be required for supporting imposed stresses. That 
is, the contact network can contain enough contacts that 
the stress balance conditions do not determine a unique 
confi gur ation of the intergrain forces. (See, for example, 
Ref. [23.) In such cases, the determination of P{f ) must 
involve some sort of average over the ensemble of pos- 
sible force configurations. Edwards has suggested that 
the appropriate measure in configuration space for this 
ensemble is a flat one; i.e., that all possible force con- 
figurations should be considered equally weighted. 23] In 
general, such an ensemble should include averages over 
different contact network geometries as well as different 



stress states on a given network. |24j 

In this paper we investigate the question of whether 
Edwards' hypothesis leads to exponential tails in P{f) 
for a system of non-cohesive grains in which the con- 
tact network forms a triangular lattice. We first study 
the case of hydrostatic compression, where our results 
confirm those of Snoeijer et al,0, 0^ though we em- 
ploy different boundary conditions and different numeri- 
cal methods. We then examine the effects of diluting the 
lattice to the rigidity percolation threshold. The transi- 
tion appears to be first order, with a finite jump in the 
number of available configurations at threshold. P(/) be- 
comes substantially broader than in the full lattice case, 
but still decays faster than exponentially on the lattice 
sizes within our numerical reach. Finally, we consider 
the effects of anisotropic imposed stresses. We show that 
strong anisotropy must produce an exponential tail in 
large systems and present numerical results showing the 
approach to this limit for varying degrees of anisotropy. 

Because the difference between the decay produced by 
our model and a true exponential decay is rather subtle, 
it is difficult to give firm interpretations of most of the 
experimental and numerical results. Observation of two 
or three decades of exponential decay could still be con- 
sistent with the slowly developing deviation we observe 
in isotropic systems. On the other hand, a given exper- 
iment may involve sufficiently strong anisotropics that 
we would expect something even closer to a true expo- 
nential. The recent results of Majmudar and Behringer 
j25l | clearly illustrate that anisotropy due to shearing pro- 
duces distributions much closer to exponential than those 
observed under isotropic compression, but we urge cau- 
tion in making direct comparisons between the model re- 
sults and experiments on disordered, deformable grains. 

In the model described below, all forces are directed 
along the line of contact between two grains and ten- 
sile forces are not allowed. This may be thought of as 
corresponding to the case of frictionless, non-cohesive. 
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circular disks, though a generic set of perfectly circular 
disks could not form as many contacts as are present in 
the model. Periodic boundary conditions are imposed, 
so there is no distinction between bulk and boundary 
contacts. By treating the force network without regard 
to any distortions of the lattice, we are considering sys- 
tems of grains with elastic moduli large enough that the 
boundary forces cause negligible strain; i.e., we consider 
either very hard grains or very weak boundary forces. In 
all cases, we fix the contact network and study the ensem- 
ble of force configurations on that network. We note that 
fiuctuations are large. P(/) in any given force configura- 
tion may look quite different from the P{f) obtained by 
averaging over configurations, even on the largest lattices 
we have studied. 

For the smallest nontrivial lattice, consisting of nine 
grains, we calculate P{f) both analytically and numer- 
ically, and find good agreement. For larger lattices, 
we employ numerical sampling methods involving Monte 
Carlo moves that maintain the force balance constraints 
at all times (in contrast to simulated annealing methods 

II. THE LATTICE 

We study the distribution of bond strengths on an 
n X n triangular lattice corresponding to the contacts 
in a hexagonal packing of monodisperse circular grains. 
We consider a system in the shape of a rhombus, having 
bonds or edges in each lattice direction, and subject 
to periodic boundary conditions, as illustrated in Fig. ^ 
Each edge carries a scalar variable / specifying the mag- 
nitude of the force transmitted across a contact. In order 
to fix the three components of the macroscopic stress ten- 
sor aij we specify the total compressive forces Fi, F2, and 
F3 supported along each of the three lattice directions. 

There are variables in the system. For each grain, 
the vector force balance condition imposes 2 constraints. 
These constraints are not all independent, however. The 
periodic boundary conditions guarantee that the sum of 
all the single-grain constraint equations is trivially zero, 
leaving 2n? — 2 independent constraints. Fixing Fi, F2, 
and imposes 3 additional constraints, for a total of 
2ti? -I- 1 constraints. This leaves — 1 degrees of freedom 
in the force configuration. 

In addition to force balance equations, there are in- 
equality constraints associated with the fact that the ma- 
terial is non-cohesive. No force / is allowed to be neg- 
ative. It is this condition that introduces nonlinearity 
in the system. For any two force balanced, tension free 



configurations on the lattice, any weighted average with 
positive weights will also be force balanced and tension 
free. Sums with a negative weighting of a configuration, 
however, will not always be allowed, as they may contain 
negative forces on some edges. 

We represent the system configuration by a vector of 
forces fi, i = 1 . . . . The above constraints create a 




FIG. 1: An nxn triangular lattice (with n = 6). Under peri- 
odic boundary conditions light edges on the left are identified 
with dangling edges on the right and light edges on top are 
identified with dangling edges on the bottom. The macro- 
scopic stress state is fixed by imposing compressive forces Fi , 
F2, and ^3 along each lattice direction. 



space of possible configurations that fills a finite, convex 
— 1-dimensional volume, with boundaries determined 
by the inequalities fi > 0. Our interpretation of the 
Edwards hypothesis is to assign a uniform probability 
(Lebesgue) measure on the space of allowed configu- 
rations. 



III. ISOTROPIC LATTICE 

In principle, P{f) can be calculated in the following 
way. Fix one edge to carry force /; this restricts the 
system to a (n^ — 2)-dimensional subset of the (n^ — 
l)-dimensional allowed volume V in configuration space. 
P{f) is simply the ratio of the (n^— 2)-dimensional "area" 
to the (n^ — l)-dimensional total volume. More formally, 

let Em be the set of 6 edges touching node to; let e^™', 
(7 = 1 ... 6, be a unit vector along edge g in Em pointing 
towards node to; and let Lfe, A: = 1 ... 3, be the set of 
edges along lattice direction k. Then we have 
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FIG. 2: A slab of material and set of edges used to prove 
that each layer of edges in Lk must have a sum of forces equal 
to Fk. 

where edge i has been fixed to have value /. The in- 
tegral is taken over all /'s. The first delta function en- 
sures that edge i carries force /. The next ensures force 
balance at each vertex, and the last enforces the bound- 
ary conditions. The integral represents the volume of a 
(n^ — 2)-dimensional slice of a polytope, therefore P{f) 
on an n X n lattice is a piecewise polynomial of order 

The average force on the edges in Lk is Fk/n. The last 
factor in Eqn. Q requires comment. For definiteness, 
consider k = 2. Let us divide the set of edges in L2 into 
n layers, each layer containing the n edges that intersect a 
line along one of the other edge directions. The following 
argument shows that the sum of the fj on each layer 
must be equal to F2 . Fig. [21 shows a shaded slab of the 
system. The total force on this slab must be zero, which 
implies that the vector sum of the forces on the edges 
indicated by thick lines on one side of the slab must be 
equal to the vector sum of the forces on the other side. 
Since the vector sum has a unique decomposition into 
contributions from the L2 edges and the L3 edges, the 
sum of the fj 's in each direction must independently be 
the same on both sides of the slab. This argument is 
independent of the thickness of the slab, so in the L2 
direction each of the n layers must have /j's that sum to 
F2, and similarly for Li and L3. 

For the 3x3 case, the integral can be evaluated analyt- 
ically. Here we state the isotropic result; the anisotropic 
result is presented below. The calculation is detailed in 
Appendix A. (We have not found an analytic expression 
for P{f) for arbitrary lattice size. See Ref. however, 
for a treatment of the case in which the sum Fi -t- F2 -I- F3 
is fixed, but not the individual Fk.) In the isotropic case 
{Fi = F2 = F3 = F), we find 

^(/)-^e(/)e(i^-/)(F-/)2 (2) 

X (5F^ -f 73/F^ - 1 11 fF^ + 125 fF^ - 59fF+9f) 

where &{f) is the Heaviside step function. This expres- 
sion is plotted in Fig. 13 

On larger lattices we employ numerical methods to 
measure P{f). This requires generating numerous con- 



FIG. 3: The wheel move associated with the node marked by 
a disk. Forces are shifted by an amount A/, chosen to be 
small enough that no edge becomes negative. Force on each 
edge marked with a thick solid segment is increased and on 
each edge marked with a thick dashed segment is decreased, 
or vice versa. 



figurations in the allowed volume of configuration space 
in a manner consistent with the Edwards fiat measure. 
Previous work has implemented a simulated annealing 
algorithm, which generates each configuration by start- 
ing from a random point in the space of all possible fi 
and relaxing to some point in the lower dimensional sub- 
set of interest. [l5lll6j | We employ a different technique in 
which, starting from one force-balanced configuration - 
a point in the allowed subset of stress states - new con- 
figurations are generated via moves that always produce 
allowed configurations. By identifying a set of moves that 
span this compact space, reach any set of positive volume 
in a finite time, and involve symmetric transition proba- 
bilities, we can be sure that at sufficiently long times the 
space is being sampled with uniform measure. |23 These 
moves are described below. 

For the nxn lattice with n> 3, we construct a set of 
wheel moves, one centered on each node. The wheel move 
associated with a given node acts on the nearest neighbor 
edges ( "spokes" ) and next nearest neighbor edges ( "rim" ) 
of the node, reducing one set of /'s and augmenting the 
other by the same amount. The move is implemented 
in the following steps. Identify the smallest force on the 
spokes and call it Smin] hkewise, call the smallest force 
on the rim Tmin- Randomly choose a force increment 
A/ with uniform measure on the interval [—Smim'Tmin]- 
Adding A/ to every edge on the spokes and subtracting it 
from every edge on the rim, which respects force balance 
on every node touching the wheel, constitutes a wheel 
move. 

In the space of linear combinations of wheel moves, 
there is one obvious null direction; i.e., a linear combina- 
tion of wheel moves that results in no change to any edge 
of the lattice. We prove in Appendix B that the number 
of linearly independent wheel moves is exactly — 1, 
which is the number required to span the space of force 
configurations. 

Fig. ^ illustrates a geometric subtlety that must be 
taken into account when attempting to explore the al- 
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FIG. 4: A two-dimensional convex space, outlined in gray, 
with two different sets of basis vectors. The basis in (a) is 
such that it is possible to move off the lower-lefthand corner 
along either basis direction, while in (b) this is not the case. 




FIG. 5: 3 X 3 triangular lattice under isotropic stress. The 
curve gives the force distribution, points represent numerical 
simulations of the same quantity. P{f) vanishes at F, the 
force imposed along each lattice direction. Here F = 3, fixing 
(/> = 1- 




FIG. 6: A typical 15 x 15 lattice. The force on an edge is 
redundantly mapped to color and width (stronger forces are 
darker and thicker). 
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lowed stress states using wheel moves (or any other basis 
set of moves) . Since the given pair of basis vectors span 
the two-dimensional convex space shown, it is clearly pos- 
sible to navigate from any interior point to any other by 
making discrete steps along the basis vectors. It may not 
be possible, however, to make a move in either basis di- 
rection if the initial configuration lies exactly on a corner. 
In higher dimensions it is likewise possible to be stuck on 
a corner or a boundary of higher dimension. It is there- 
fore important to find initial configurations that lie in 
the interior of the space, rather than on a boundary, for 
the purposes of our Monte Carlo sampling technique. An 
algorithm such as the one described above, beginning at 
an interior point, can come arbitrarily close to all bound- 
aries but will never reach them. Regions near corners are 
visited infrequently but for long times in such a way that, 
for sufficiently long runtime, the space is sampled with 
uniform measure. 27] 

Fig. 13 shows the agreement between the exact calcu- 
lation and simulation on the 3x3 lattice. The forces are 
normalized by choosing F such that the average force (/) 
is unity. The peak near (/) is typical for larger lattices 
as well, but finite size effects are clearly evident in the 
tail, since P{f) must go to zero at f = F (here F = 3). 

A typical configuration for the 15x15 lattice is shown 
in Fig.|Sl -P(/)'s for the cases n — 5, 10, 15, and 20 are 
shown in Fig.Q The peak near (/) is again apparent, and 
for small / the curves appear to coincide, though small 
differences exist that are masked by the logarithmic scale. 



FIG. 7: Force distributions for n x n lattices under isotropic 
stress for n = 5,10,15, and 20. For all cases the decay is faster 
than exponential. The n = 15 and n = 20 cases are nearly 
identical, indicating convergence to a universal curve for large 
systems. The straight line is drawn as a guide to the eye. For 
each curve, (/) = 1. 

For / 3(/) the curves separate, with the n = 5 distri- 
bution decaying most rapidly. The four distributions in 
the figure all decay faster than exponentially, and they 
broaden slightly with increasing n. 

To gain confidence that the curves are converging to 
a large system limit, we measured correlations between 
forces as a function of the distance between edges. In di- 
rections both longitudinal and transverse to a particular 
edge, as shown in Fig.|S| we see that on a 20 x 20 lattice 
correlations have decayed to the 1% level at a distance of 
10 lattice constants. On larger lattices we see exponential 
decay with a decay length of approximately 3 lattice con- 
stants, although much more data on larger lattices would 
be necessary to measure this precisely. We therefore ex- 
pect to see little difference in P{f) for lattices larger than 
n — 20. Supporting this expectation, there is very little 
difference between the n = 15 and n = 20 distributions 
for / up to the largest values we have measured, which 
covers five decades of P. Though it is conceivable that 
the asymptotic form of P{f) at large / in our model is 
exponential, the domain displayed in Fig. is the rele- 
vant one for comparison with experiments, and it clearly 
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FIG. 8: Above, correlations are calculated in the longitu- 
dinal and transverse directions. Below, correlation function 
g{r) = (fifi^r) — (fi) on a 20 X 20 lattice. The correlation 
is computed for edges that have the same orientation. The 
longitudinal correlation refers to displacements by r lattice 
constants along the direction of the edge. The transverse 
correlation refers to displacements perpendicular to the edge. 
The plots above are averaged over multiple reference edges. 
The transverse correlations are in fact negative, but we take 
their absolute value to facilitate plotting on a log scale. 



shows a decay that is faster than exponential. 



IV. DILUTED LATTICES 

Real disks are not perfectly monodisperse, so in any 
hexagonal packing of hard disks some of the bonds on 
the triangular lattice will not actually be present. For 
a given stress state, a diluted lattice has fewer bonds to 
carry the same force. Those remaining will, on average, 
carry more force, resulting in a shift in P(/) toward larger 
forces, which may involve a broadening of the tail. As the 
number of edges removed from the lattice increases, the 
force configurations tend to become less homogeneous, 
with strong forces concentrated in chain-like structures. 
As more edges are deleted there comes a point where the 
lattice can no longer support the imposed stress; this is 
the rigidity-percolation transition. Fig.|51shows a typical 
force configuration for a randomly diluted lattice. We 
study randomly diluted lattices at the transition to see 
whether the broadening takes the form of an exponential 
decay. 

The random dilution process merits further remark. 
For a random process in which each bond is removed 




FIG. 9: A 15x 15 lattice at rigidity percolation. Deleted edges 
are covered with disks. The force on an edge is redundantly 
mapped to color and width (stronger forces are darker and 
thicker) . 



with probability cf), the infinite triangular lattice un- 
der isotropic compression cannot support stress for any 
(j> > 0. j23| Thus in an rt X n simulation the fraction 
(not the number nd) of deleted edges at the rigidity per- 
colation threshold goes to zero as n — > oo. As an arbi- 
trarily large real hexagonal packing does support com- 
pression, this means that the process of bond breaking 
in the real packing does not happen randomly but rather 
in a correlated way. Nevertheless, our simulations em- 
ploy random dilution on finite lattices. We also neglect 
for present purposes the possible buckling instability of a 
force chain. Disallowing configurations that might buckle 
could only result in ensembles with fewer dilutions and 
therefore would not change the conclusions described be- 
low. 

We construct lattices at the rigidity percolation thresh- 
old in the following way. Given the imposed macroscopic 
stress and beginning with an empty lattice, edges are 
selected at random to be added to the lattice. After 
each addition, we check to see if the lattice is capa- 
ble of supporting the imposed stress using the simplex 
algorithm. 29] We take the lattice to be at threshold when 
it first supports the imposed stress, at which point we 
know that there exists at least one edge such that, when 
removed, the system could no longer support the stress. 
The edges that have not been added to the lattice at this 
stage are called "deleted edges." Our construction yields 
a lattice in which deletions are randomly distributed in 
space. Note that the process will in general yield some 
edges that are present but not connected to the edge net- 
work in a way that allows them to bear any force. We 
refer to these edges as "effectively deleted." 

These effectively deleted edges on a lattice supporting 
compressive forces can be identified in the following way. 
Consider the 6 edges that meet at a given node. An edge 
becomes effectively deleted if its opposite edge has been 
deleted (possibly effectively) and at least one of the edges 
making a 120° angle with it has also been deleted (possi- 
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bly effectively). Under these conditions any force on the 
edge in question can not be balanced by positive forces 
on the other edges sharing the node. To find all of the 
effectively deleted edges we examine all nodes repeatedly 
until no new deletions are found. 

The wheel moves used for investigation of the undi- 
luted lattice no longer work on the diluted lattice. The 
wheel moves add or subtract a quantity A/ to each edge 
they touch, but this is impossible if one of the edges is 
deleted. Therefore each deleted edge renders unusable 
the four wheel moves to which it belongs. To salvage a 
set of moves that span the appropriate space, linear com- 
binations of the wheel moves can be formed that leave 
the deleted edge unaffected. For a single deleted edge 
there are three linearly independent combinations of the 
four wheel moves that edge touches which preserve the 
deleted edge. The deletion has reduced the available de- 
grees of freedom by one. When many edges have been 
deleted, more complicated linear combinations of wheel 
moves can be found that preserve the vanishing force on 
all deleted and effectively deleted edges. 

We refer to a particular linear combination of wheel 
moves that remains a viable move on a diluted lattice as 
a "multi-wheel move." Let the number of linearly inde- 
pendent multi- wheel moves be Nm- Just as the undiluted 
lattice had in? wheel moves and n? — \ degrees of free- 
dom, the diluted lattice having Nm multi-wheel moves 
has Nm — 1 degrees of freedom. Finding a complete set 
of multi-wheel moves then provides not only a means to 
perform numerical simulation but also a count of the de- 
grees of freedom in the diluted system. 

A linearly independent set of Nm multi-wheel moves 
on the diluted lattice can be constructed as follows. For 
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form a vector Hi j £ R"'' such that the 



components of l6 j specify the effect of the wheel move 
centered on node j on the deleted edges. That is, {v^j)i — 
+1 if the i*"^ deleted edge is a spoke of node j, —1 if 
the i"^ deleted edge is a rim of node j, and otherwise. 
Form the Ud x n? matrix Wij — Consider a linear 

combination of wheel moves with coefficients m^, k = 
1 . . .n^. This linear combination has no effect on deleted 
edges if and only if Wfft = 0. Thus 7V,„ = dimiV(Vl^) = 
r? — rank(iy), where N{W) is the null space of W. 

Note that we demand the wheel moves respect both 
real and effective deletions. Were we not to specify effec- 
tive deletions as well, the algorithm would in general pro- 
duce additional multi-wheel moves which were nonzero 
on some effectively deleted edges. This is because the 
multi-wheel moves could also act on a lattice that sup- 
ports tensile as well as compressive forces. On such a 
lattice there are fewer effective deletions; in particular, a 
node could be force-balanced having only three edges all 
on the same side of a line through the node. To model a 
noncohesive material, then, we must include these addi- 
tional effective deletions as well. 

We use the following technique to investigate the de- 
grees of freedom in lattices at threshold. Edges are 
added to a lattice one at a time until the rigidity per- 
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FIG. 10: The ratio of average number of degrees of free- 
dom in threshold lattices to the number in the undiluted case 
(S) = ((iVm) - - 1) vs. the ratio of deleted ed ees in the 

threshold lattice to the total number of edges (j) = Udj'ir^ ■ 
For each plot a number of lattices at the rigidity-percolation 
threshold were generated and the degrees of freedom and 
deleted edges counted. For large enough lattices, the system 
has a large number of degrees of freedom available as soon as 
rigidity-percolation is reached. Note that the observed inter- 
val of A values shifts as the lattice size is increased. 



eolation threshold is reached, the set of wheel moves 
on that threshold lattice are constructed to determine 
the dimension of the volume of allowed configurations 
in state space Nm — 1, and the process is repeated for 
a number of threshold lattices. We refer to the sys- 
tem as having Nm — 1 degrees of freedom, as this is the 
number of multi-wheel moves it possesses. These lat- 
tices have differing numbers of deleted edges. For lattice 
with the same number of deleted edges rid we calculate 
the average ratio of the number of degrees of freedom 
of the dilute system to the number in an undiluted one: 
(8) = {{Nm) - - 1). Figure Cni shows (5) plotted 

against the fraction of of deleted edges <^ = ndjiv? . A 
discrete jump in {S) emerges as is decreased for sys- 
tems larger than 15 x 15, indicating that the transition is 
first order. This suggests that near the transition 

point is unlikely to show qualitatively different behavior 
from that of the undiluted lattice, though a quantitative 
broadening of the force distribution is expected. 

We measure -P(/) averaged over a number of thresh- 
old 20 X 20 lattices. The resulting distribution is indeed 
much broader than the undiluted case, but there is still 
curvature in the distribution on a linear- log plot. Numer- 
ical studies in this regime are hampered by two factors. 
First, each multi- wheel move requires the examination of 
many edges to ensure that no /j becomes negative. Sec- 
ond, the maximum size of a multi-wheel move is typically 
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FIG. 11: P{f) for 20 x 20 triangular lattices diluted to the 
rigidity-percolation threshold. The plus signs represent an 
average over 2.5 x lO'^ moves on a single randomly diluted 
lattice. The open circles represent an average over 75 lattices, 
each run for 5 x 10* moves. The straight line is a guide to the 
eye. 



quite small because the move touches several edges that 
carry small forces. The volume that we are attempting to 
sample with uniform measure is still convex but contains 
many corners that are not easily accessed. The data ob- 
tained for a given lattice may exhibit a bias depending 
on the initial configuration unless the number of moves 
considered is extremely large. 

We have gathered data using two procedures. In one 
case we construct a large number of threshold lattices, 
find a single initial configuration for each and collect force 
data from 5 x 10^ wheel moves, then average all of the 
data together to determine P{f) averaged over lattices. 
To find initial configurations we use the simplex method 
several times with different coefficients and average the 
results, thereby avoiding configurations that lie on the 
boundary of the allowed volume. We cannot be sure, 
however, that this method is free of systematic bias. In 
the second case, we consider a single lattice at threshold 
(with a typical number of degrees of freedom), find 25 dif- 
ferent initial configurations and run 2.5 x 10^ multi-wheel 
moves from each one, averaging the force data to obtain 
the P{f) associated with that particular lattice. The re- 
sults are shown in Fig. ^] As may be expected given the 
first-order nature of the rigidity percolation transition, 
neither procedure produces an exponential tail in P at 
large /. We conclude that random dilution to the rigidity 
percolation threshold does not produce exponential tails 
for triangular lattices. 



V. ANISOTROPIC LATTICE 

We now consider the undiluted lattice with anisotropic 
stresses imposed by choosing Fi different from F2 and 
^3. For example, one lattice direction may be subject to 
stronger compression than the others, creating qualita- 
tively distinct force distributions in the strong and weak 
directions. We will show that in strongly anisotropic sys- 




FIG. 12: a — > 00 limit of the nxn triangular lattice. Dashed 
edges represent attachments under periodic boundary condi- 
tions. The macroscopic stress state is fixed by imposing com- 
pressive force A along the strong lattice direction, the only 
one carrying force. Force is transmitted along n chains, each 
composed of n edges, in any way such that a sum over chains 
yields A and all forces remain compressive. 



tems the strong direction contributes an exponential de- 
cay to the tail of P(/). 

We consider Fi = F + A, F2 = = F, and pa- 
rameterize the anisotropy by a = (F -|- A)/F, so that 
a = 1 corresponds to isotropic stress. In terms of (/), 
F = 3n(/)/(2 -I- a). The average force in the weak di- 
rection {fm) = F/n, while that in the strong direction is 
(/s) = A)/n. Similarly, forces in the weak direction 
can not exceed F, and those in the strong direction can 
not exceed F -I- A. Forces in the strong direction will 
populate the tail of P{f)- 

In the limit a — > 00 the system becomes much simpler: 
only the strong direction carries force, doing so via n 
chains of n edges, each edge in a given chain carrying 
the same force as shown in Fig. ^| The forces on the n 
chains must sum to A. 

Simple dimensional analysis of the scaling of the al- 
lowed volume when one chain is fixed at force / yields 



lim P{f) oc lim 1 — 



/ 



n{fs 



-f/ifs) 



(3) 



We see that in this limit anisotropy induces an exponen- 
tial tail in P{f)- We next investigate the extent to which 
the limiting behavior is reflected at finite anisotropies. 

To gain some intuition for the approach to the 
anisotropic limit we return to the 3x3 case, which can be 
solved analytically. We let Ps{f) and Pw{f) denote the 
separate distributions of forces in the strong and weak 
directions, respectively. Pw{f) is identical to P{f) in 
Eqn. 13 In the strong direction, however, the situation 
is more complicated. A complete expression is given in 
Appendix C. 

For the 3x3 case, the fully anisotropic limit a — > 00 
has P[f) = |(/,)-2(3(/,) -/), / G [0,3(/,), a piecewise 

linear function. For a > 2, P{f) is given by Pi^^\f) 
(from Appendix C) between F and A. pi^^^ is a linear 
function, and the interval from F to A grows with a. As 
a increases further the width of this region grows and the 
slope approaches the limiting slope of — |(/s)^^. Ps{f) 
is plotted for increasing a in Fig. 1131 



8 




FIG. 13: 3x3 triangular lattice under anisotropic stress. The 
force distribution Ps{f) for forces in the strong direction is 
pfotted for increasing anisotropy. Each distribution is scaled 
to (fs) = 1. The limiting form is a piecewise linear function 
of /. Pa for finite anisotropy has a linear region in the middle; 
this region grows in width and approaches the limiting slope 
as anisotropy is increased. 




FIG. 14: A 15 X 15 anisotropic lattice with a = 5. The 
force on an edge is redundantly mapped to color and width 
(stronger forces are darker and thicker). The stronger stress 
is imposed in the horizontal direction. 



On larger lattices, numerical calculations show similar 
behavior to the 3x3 case, but with the limiting linear 
distribution replaced by the limiting exponential of the 
large system limit over the range of forces of interest. 
(The distribution must be cut off at / = Fg.) Fig. H"!] 
shows a typical anisotropic configuration. Fig. 1151 shows 
P{f) for anisotropic lattice of size 15 x 15, with the con- 
tributions from the weak and strong directions shown 
separately. The role of the strong direction in broaden- 
ing the distribution is clear. Results are shown for two 
cases, one in which strong forces are imposed on a single 
lattice direction and another in which strong forces are 
imposed on two lattice directions. Note the difference in 
behavior for small /. In the latter case there is no peak 
in P(f) for small /. 

Fig. 1 161 shows force distributions on a 15 x 15 lattice for 
several values of a. For sufficiently large a a middle por- 
tion of P{f) appears to be exponential. As a increases, 
so does the extent of this portion. We hold (fs) fixed so 



A weak (2) 
+ strong 
all 




FIG. 15: Pi.f) for strong, weak, and all lattice directions 
on a 15 X 15 lattice with a = 3. The strong direction(s) is 
entirely responsible for the tail in P{f)- Top: Two weak J^'s 
and one strong. Bottom: One weak Fk and two strong. On 
both plots, (/) = 1. 
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FIG. 16: P{f) on a 15 X 15 lattice for increasing anisotropy. 
Portions of the tail for stronger anisotropies are difficult to 
distinguish from exponential decay, and the limiting behavior 
is a true exponential (gray line). The dashed line indicates 
the interval in / over which the a = 8 curve is approximately 
exponential, (/s) = 1 for every curve. 



that every curve will have the same decay length. 

We conclude that P{f) on the anisotropic lattice dis- 
plays an exponential tail in the following sense. For mod- 
erate a and n, a portion of the tail of P{f) is nearly linear 
on a log plot. The decay is not truly exponential because 
of the finite anisotropy and finite lattice size, but the 
limiting behavior of P{f) as n and a arc increased is a 
true exponential. The region of nearly exponential de- 
cay grows with increasing anisotropy, extending as far as 
/ « 3(/,) for a = 8. 

We note that the degree of anisotropy represented by 
a given value of a can be compared to the anistropy sup- 
ported by materials described by an internal friction pa- 
rameter. The ratio of major and minor principal stresses 
on our lattice is cri/a'2=(a -|- 1)/^. By exploiting the 
relation axjai = (1-1- sin (/))/( 1 — sin 0) [S^l , where 4> is the 
angle of internal friction, we find that a = 8 corresponds 
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to an internal friction of approximately 43°. 

The survival of the exponential tail of strong forces 
in large triangular lattices under anisotropic loading can 
be traced to the discussion based on Fig. |21 concerning 
the sums of forces along layers of edges in the various 
directions. That argument shows that the strong forces 
in one lattice direction cannot be redistributed into the 
other directions due to vector force balance constraints. 
Thus for strongly anisotropic lattices, the strong forces 
will follow scaling laws close to the limiting exponential 
form corresponding to no force at all in the weak direc- 
tion(s). Note that in the case of two strong directions 
and one weak, the two strong are effectively independent 
since strong chains from the two directions cannot in- 
teract significantly without generating a strong chain in 
the supposedly weak direction and thereby reducing the 
degree of anisotropy. 



VI. CONCLUSION 

We have investigated the distribution P{f) of contact 
forces on a lattice of triangular bonds under the Edwards 
flat measure. An interesting question is whether one ex- 
pects exponential decay in P{f). The triangular lattice 
is a simple but nontrivial system for studying this phe- 
nomenon. 

The distribution on the n x n periodic lattice decays 
faster than exponentially, as reported previously. ^ Di- 
luting the lattice induces significant broadening in P(/), 
but the decay remains faster than exponential. Even at 
rigidity percolation associated with random bond dilu- 
tion of the lattice, no qualitative change in the form of 
P{f) is discernible. In particular, we do not see evi- 
dence for an exponential tail associated with the transi- 
tion. Consistent with these direct measurements of P{f), 
we find that the transition is first order, which generally 
suggests that no qualitative changes should be expected 
as the percolation threshold is approached. 

On the other hand, imposing anisotropic stress on the 
undiluted lattice can induce an exponential tail in P{f). 
In the limit of an infinite lattice with stress imposed only 
along one lattice direction, the distribution of contact 
forces is a pure exponential. Numerical simulation shows 
that evidence of this behavior may still be seen for finite 
lattice sizes and a finite ratio of the compressive forces in 
the strong and weak directions. In such a scenario the tail 
of P{f) is not a true exponential, but appears approxi- 
mately linear on a log plot of P{f) for some interval in /. 
For large enough anisotropies, P{f) decays three orders 
of magnitude from its maximum and / w 3(/s) before 
deviation from exponential decay becomes obvious. 

In the triangular lattice model, anisotropy is associated 
with the appearance of long force chains oriented along 
the strong direction, a phenomenon that has also been 
observed in experiments on disordered svstems.|25l iBlj 
Extension of the numerical techniques employed above 
to rigid bars that form a disordered triangulation of the 




(b) (c) 

FIG. 17: Three of the nine basis elements employed to evalu- 
ate Eqn.Qfor the 3x3 isotropic lattice. Solid bars indicate 
a positive contribution to the force on an edge; dashed bars 
indicate a negative contribution. Double bars indicate a con- 
tribution with twice the weight. The honeycomb in (a) is one 
of three. The elements in (b) and (c) are each repeated for 
the other two lattice directions. 



plane would be straightforward. Wheel moves, for exam- 
ple, would be defined by taking all bars sharing a given 
vertex to be the spokes, though it would no longer be 
true that all bars receive force increments of the same 
magnitude during a move. The primary difference be- 
tween the prefect triangular lattice and disordered ones 
(or even other crystalline ones) is that lattices with no 
sets of collinear bars that span the system cannot sup- 
port arbitrarily strong anisotropic stresses; i.e., in most 
cases there would be a maximal value of a. A detailed 
investigation of the effects of disorder would be of inter- 
est. 
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APPENDIX A: ANALYTIC P{f) 

P{f) for the 3x3 isotropic lattice is calculated as 
follows. There are 3^ — 1 = 8 degrees of freedom; we 
take nine basis elements subject to a constraint which 
is to be imposed by hand. Three of these elements are 
shown in Fig. 1171 There are three honeycomb elements, 
i = 1 . . .3, three elements i j = 1 • ■ ■ 3, and three 
elements ■0^2 • The figure depicts i/'ii, and V'i2- The 
ip elements make no net contribution to the total force in 
the system. The three 4> elements must sum to F; this 
is the additional constraint we impose. All elements are 
isotropic by construction. 
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FIG. 18: Useful references for the proof in Appendix B. 



Eqn.^can be rewritten 



P{f) 



N 



n 



bo-4']2 



X <5(/ - (01 - 21/^12)) 

X S{(Pi + (l)2 + (l)3- F) 



(Al) 

where (f)o = min({(/)i, 02, ^a})- The first (5-function fixes 
a particular edge to support force /; all edges are equiva- 
lent by symmetry. The bounds on the ip^s assure that no 
edge will support tensile force. The presence of 0o breaks 
the integral into three regions. The regions wherein 
00 — 02 and 00 = 03 are identical; the region where 
00 — 01 must be evaluated separately because 0i touches 
the edge fixed to carry /. 

The evaluated expression is given in Eqn. |31 



APPENDIX B: PROOF THERE ARE 1 
INDEPENDENT WHEEL MOVES 

Adapting our notation from Section IV, for each j — 
\ . . .n? let j e R"^" be the vector that specifies the 
effect of the wheel move centered at node j on the collec- 
tion of all edges. Form the iv? x v? matrix Wtj — {nS 
A linear combination of wheel moves with coefficients 
mfc, k = 1 . . .v? has no effect on any edge if and only if 
Wrft = 0. 

Lemma: N[W) is spanned by a single vector such that 
rrii = rrij for all i, j. 

Remark: Since W is x n^, the system of equations 
Wm = is highly overdetermined. 

Proof: Given any edge, let i and j be the two vertices 
at the ends of the edge, and let k and I be the nearest 
two vertices on the perpendicular bisector (see Fig. 118b). 



The wheel move with coefficients {mi\ has no effect on 
the given node if and only if 



rrii + rrij — rrik — mi = 0. 



(Bl) 



To interpret Eon. IBll geometricallv. let {^i,Vi) be the 
coordinates of node i in Fig. 118b. and define Xi £ R'' by 

= (/ij,z^j,mj). (B2) 
Define Xj , Xk , and Xi similarly. Then Eqn. IBll is satis- 
fied if and only if Xi lies in the plane in R'' determined 
by X^, Xj, and Xk- 

Suppose the wheel move with coefficients {rrii} has no 
effect on any edge: i.e., Wrft = 0. Order the nodes 
as indicated in Fig. 118b . For the moment we regard the 
first three coefficients mi , m2 , as arbitrary. Applying 
Eqn. IBll to nodes 4, 5, and 6, we deduce that X4, X^, 
and Xq all lie in the plane of Xi, X2, and X3. 

With further application of Ean. lBll we can extend the 
conclusion to all vertices of the outer triangle, and this 
may be continued to conclude that all points Xi lie in a 
single plane. Finally, periodicity requires that this plane 
is horizontal: i.e., rrii = for all i and j. QED 



APPENDIX C: ANISOTROPIC P{f) ON THE 3x3 
LATTICE 

The calculation of P{f) for the anisotropic 3x3 
case proceeds similarly to the isotropic calculation. The 
anisotropy requires the division of the relevant integral 
according to whether the force / is greater than or less 
than the difference A between the strong Fi and the weak 
ones. For the case A < F, we find 



P?"V) A<f<F 
Pf^(/) F</<F + A, 



Psif) 



1 

N 



(CI) 



where 1^ is a normalization constant and the Ps's are 
given below. For the case A > F 



p^'Hf) 

p^''\f) 
5(3) 



Psif) = ^ 



The functions Ps in Fan. IC2l are 



f<F 
F<f<A 



(C2) 



Pr(/) A</<F + A. 



pji)(/) = (7a-2)P^ + 21(4a-l)PV-42(5a-|-l)PV^ + 140(2a + l)P^/^ 

-210(a + l)P^f'* + 84(a + 2)F'^f^ - 14(a + b)F f^ + 12/^; (C3) 
p(^^)(f) = (3a^ - 14a^ + 21a^ - 35a^ + 42a^ - Ua + 2)P^ -I- 21 (a'^ - 4a^ -1- 5a^ - h)a^F^f 
+2l{?,a^ - lOa^ + lOa^ - 15)aPV^ - 35(3a^ - Sa^ -I- Qa^ -8a- 5)PV^ 



11 



+105(^3 -2a^ -a- 2)F^f - 21{3a^ -8a- 7)F'^f + 7{a - 12)Ff + 9f; (C4) 

^ (21a - 4)F^ 21F^/; (C5) 
Pj3)(/) = (aF - ff [(-3a^ - 14a'* - 21a=^ + 35a - 42)F^ + (ISa'' - SGa^ + GSa^ - 35)FV 

+(30a2 - 56a + 21)i^V^ - 3(100^ - 28a + 21)af - (15a - 14)^^/-* + 3/^] . (C6) 
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